
% Position specific frequency matrix (PSFM) for TBP/TATA-box from JASPAR.
M1 = [61   16  352    3  354  268  360  222  155   56   83   82   82   68   77;
     145   46    0   10    0    0    3    2   44  135  147  127  118  107  101;
     152   18    2    2    5    0   10   44  157  150  128  128  128  139  140;
      31  309   35  374   30  121    6  121   33   48   31   52   61   75   71];

M1 = M1 + ones(size(M1)); % Add a pseudo count.
M1 = M1./(ones(4,1)*sum(M1,1)); % normalize.
M_1 = {M1};

% PSFM for AP2 (TFAP2A) from JASPAR.
M2 = [  0   0   0  22  19  55  53  19   9;
        0 185 185  71  57  44  30  16  78;
      185   0   0  46  61  67  91 137  79;
        0   0   0  46  48  19  11  13  19];
M2 = M2 + ones(size(M2)); % Add a pseudo count.
M2 = M2./(ones(4,1)*sum(M2,1)); % normalize.
M_2 = {M2};
M_12 = {M1,M2};

% NFYA / CAAT-BOX
M3 = [34  16   7  58  51   0   2 112 116   0  14  66  13  39  36  25;
      37  33  51  14   4 116 113   0   0   1  65   6  20  43   9  35;
      27  26  25  41  56   0   1   1   0   0  33  42  73  22  47  29;
      18  41  33   3   5   0   0   3   0 115   4   2  10  12  24  27];
M3 = M3 + ones(size(M3)); % Add a pseudo count.
M3 = M3./(ones(4,1)*sum(M3,1)); % normalize.
M = {M3};

%SP1 / ZN-FINGER
M4 = [1 2 0 0 0 2 0 0 1 2;
      1 1 0 0 5 0 1 0 1 0;
      4 4 8 8 2 4 5 6 6 0;
      2 1 0 0 1 2 2 2 0 6];
M4 = M4 + ones(size(M4)); % Add a pseudo count.
M4 = M4./(ones(4,1)*sum(M4,1)); % normalize.
M = {M4};

% A Markovian background model: zero:th order
Bd = [0.2799 0.2201 0.2201 0.2799];

% Read sequence
[S,L] = readfastaseqs('Myod1.fa');
S = S{1};
S = basepairs2num(S);

% Read additional data: conservation scores.
[PC,L] = readPhastConsScoresSingle('Myod1-Merged.gff');
PC = PC{1};
PC = (2*0.1)*PC + 0.5 - 0.1; % Soften the prior/scale between 0.4 and 0.6.

% Run: TATA-box
[Pc1,c_mu1,priorA1] = TFBS_DF_Likelihood(S,[],M_1,Bd,0);
p1 = 1 - Pc1(1)
Pc1(1:5)
[Pc2,c_mu2,priorA2] = TFBS_DF_Likelihood(S,PC,M_1,Bd,0);
p2 = 1 - Pc2(1)
Pc2(1:5)

% Run: AP2
[Pc3,c_mu3,priorA3] = TFBS_DF_Likelihood(S,[],M_2,Bd,0);
p3 = 1 - Pc3(1)
Pc3(1:5)
[Pc4,c_mu4,priorA4] = TFBS_DF_Likelihood(S,PC,M_2,Bd,0);
p4 = 1 - Pc4(1)
Pc4(1:5)
pp1 = p1*p3
pp2 = p2*p4

%return;

% Binding site prediction at single nucleotide resolution.
% NOTE: input 'priorA' comes from 'TFBS_DF_Likelihood'.
[postnm1,postloc1,postcomb1,totalsamples1] = ...
    TFBS_DF_Bayes_MCMC_convdiag(S,[],0,Bd,M_1,1,priorA1,[],...
    100000,100000,1000000,0.05,2);
figure;
subplot(211);
area([1771 1777],[1 1])
hold on
plot(postloc1)
axis([1 2000 0 1])
ylabel('binding prob.')
xlabel('position relative to TSS');
%figure;
subplot(212)
area([1771 1777],[1 1])
hold on
plot([1751:1:1800],postloc1(1751:1:1800))
axis([1751 1800 0 1.01])
ylabel('binding prob.')
xlabel('position relative to TSS');

[postnm2,postloc2,postcomb2,totalsamples2] = ...
    TFBS_DF_Bayes_MCMC_convdiag(S,PC,0,Bd,M_1,1,priorA1,[],...
    100000,100000,1000000,0.05,2);
figure;
subplot(211);
area([1771 1777],[1 1])
hold on
plot(postloc2)
axis([1 2000 0 1])
ylabel('binding prob.')
xlabel('position relative to TSS');
%figure;
subplot(212)
area([1771 1777],[1 1])
hold on
plot([1751:1:1800],postloc2(1751:1:1800))
axis([1751 1800 0 1.01])
ylabel('binding prob.')
xlabel('position relative to TSS');


[Pc,c_mu,priorA] = TFBS_DF_Likelihood(S,[],M_12,Bd,0);
[postnm3,postloc3,postcomb3,totalsamples3] = ...
    TFBS_DF_Bayes_MCMC_convdiag(S,[],0,Bd,M_12,1,priorA,[],...
    100000,100000,1000000,0.05,2);
[postnm4,postloc4,postcomb4,totalsamples4] = ...
    TFBS_DF_Bayes_MCMC_convdiag(S,PC,0,Bd,M_12,1,priorA,[],...
    100000,100000,1000000,0.05,2);

% Analyze both strands.
[Pc5,c_mu5,priorA5] = TFBS_DF_Likelihood_Double(S,[],M_1,Bd,0);
1 - Pc5(1)
Pc5(1:5)
[Pc6,c_mu6,priorA6] = TFBS_DF_Likelihood_Double(S,PC,M_1,Bd,0);
1 - Pc6(1)
Pc6(1:5)

